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We present results from our study of the Parallel Tempering algorithm. We examine the swapping acceptance 
rate of a twin subensemble PT system. We use action matching technology in an attempt to maximise the 
swap acceptance rate. We model the autocorrelation times within Parallel Tempering ensembles in terms of 
autocorrelation times from Hybrid Monte Carlo. We present estimates for the autocorrelation times of the 
plaquette operator. 


1. INTRODUCTION 

Dynamical fermion simulations are still very 
demanding computationally. While a modern su¬ 
percomputer is able to produce 0(100) quenched 
gauge configurations in an afternoon, months are 
needed to generate a similar number of dynamical 
configurations on the same machine. 

Tempering algorithms have been successful in 
the past in reducing autocorrelation times in dif¬ 
ficult spin glass simulations such as the random 
field Ising model [|j. It was of interest to see 
whether parallel tempering (PT) could achieve 
the same success for lattice QCD. 

We present only our main results in this pub¬ 
lication. Full details of our work can be found in 
Hj. A lot of the background information is cov¬ 
ered in [||-||. 

In section ^ we outline the PT algorithm. In 
section ^ we present the predictions of a simple 
model of autocorrelations in a twin subensemble 
PT system. Our simulation parameters are out¬ 
lined in Results from the simulations are pre¬ 
sented in section ||. Our summary and conclu¬ 
sions are in section || 

2. THE ALGORITHM 

Parallel tempering consists of running several 
simulations, hereafter referred to as subensem¬ 
bles, concurrently. Each subensemble has its own 
Hamiltonian T-Li, parameter set and phase space. 
The overall PT state then is the set {si|i = 1...A} 
where Si is the state of subensemble i and N is the 


number of subensembles. The PT phase space is 
the direct product of the individual phase spaces. 

The goal of PT is to construct a Markov Pro¬ 
cess which will converge to the equilibrium prob¬ 
ability distribution 

( 1 ) 

I 

where the terms on the right hand side are 
the equilibrium distributions of the individual 
subensembles at their given parameter sets. 

One defines two kinds of Markov transitions: 

• transitions within subensembles 

• transitions between subensembles 

Transitions within subensembles are carried 
out using HMC. Transitions between subensem¬ 
bles involve a proposal to swap the current fields 
in a subensemble i with those in subensemble j. 
The swap proposal is accepted with the Metropo¬ 
lis acceptance probability 

Psii,j)=nim{l,e~^'^), (2) 

where 

An = {nj(a) + n^(b)} - {n^(a) + n^m , (3) 

which satisfies detailed balance with respect 
to Peg by construction. The resulting overall 
Markov Process is connected and satisfies detailed 
balance with respect to the required equilibrium. 


2.1. Swap Acceptance 

The acceptance rate of swap proposals deter¬ 
mines any reduction in autocorrelation times over 
the usual HMC ones. It is also related to a 
distance in parameter space via the formalism 
of action matching technology |^. It can be 
shown from detailed balance considerations that 
the swap acceptance rate is given by 

(A)=erfc(ix/^). (4) 

3. AUTOCORRELATIONS 

Consider a twin subensemble system where 
both subensembles have the same HMC auto¬ 
correlation function Cuit). This assumption is 
justifiable if the HMC autocorrelation times vary 
slowly over the region of parameter space where 
one intends to carry out PT simulations. Table 
1^ shows the integrated autocorrelation times (th) 
of the plaquette for two HMC simulations sepa¬ 
rated by such a distance in parameter space. It 
can be seen that the integrated autocorrelation 
times are equal within errors. 

Let us assume that individual swap probabili¬ 
ties may be replaced by the average swap prob¬ 
ability (A). Furthermore we are interested only 
in an even number of swap attempts. Only after 
an even number of successful swaps can a config¬ 
uration end up in its original subensemble. We 
assume no cross correlations between subensem¬ 
bles. 

It can be shown Q that the resulting PT au¬ 
tocorrelation function in each subensemble is: 

CpT{t) = \{l + {l-2{A)f]cH{t) (5) 

and assuming an exponentially decaying HMC 
autocorrelation function, the ratio of the HMC 
integrated autocorrelation time to its PT coun¬ 
terpart (tpt) is given by 

TPT ^ {A)(th - 1) -fl , , 

Th 1 + 2{A)th 

The above ratio is bounded above by 1 and tends 
to i as Th increases for a fixed (A). Hence this 
model predicts that the expected gain from PT 
over HMC is a factor of 2 in each subensemble 


Table 1 


HMC integrated autocorrelation times for the 
plaquette 


/3 

K, 

Th 

5.2 

.1335 

18(8) 

5.232 

.1335 

20(6) 


for roughly twice the work. Models of systems 
with different HMC autocorrelation times are un¬ 
der investigation. 

4. SIMULATION 

Our simulations used the GHMC code devel¬ 
oped by the UKQCD collaboration for HMC tran¬ 
sitions with extra logic to carry out the swapping. 
The simulation parameters to be tuned were the 
inverse gauge coupling (3, the fermion hopping pa¬ 
rameter K and the clover coefficient c. 

We carried out five twin subensemble simula¬ 
tions: 51, 5'2, 53, S'4 and 5'5. The first subensem¬ 
bles of each one had parameters (/3i,ci,«:i) = 
(5.2,2.0171,0.1330). The parameters for the sec¬ 
ond subensembles are shown in table ||. The sim¬ 
ulation parameters for 51, 52 and 53 were tuned 
using action matching [|j. For 54 and 55 only k 
was varied. The lattice size used was 8^ x 16. A 
reference HMC run at the parameters of the first 
ensembles was also carried out. 

5. RESULTS 

Table || summarises the results of our simu¬ 
lations. Looking at columns 2, 3 and 4 it can 
be seen that (AH) is usually greater than one 
and that (A) drops rapidly as the magnitude of 


Table 2 

Simulation parameters for second ensembles 


Simulation 

(/32,C2,K2) 

51 

(5.2060,2.01002,0.13280) 

52 

(5.2105,2.00471,0.13265) 

53 

(5.2150,1.99940,0.13250) 

54 

(5.2,2.0171,0.13280) 

55 

(5.2,2.0171,0.13265) 











Table 3 

Simulation Results 


Simulation 

Ak(x10-^) 

(AH) 

(^) 

Ant 

TH 

HMC 

- 

- 

- 

26(6) 

1 

51 

-2.0 

1.23(2) 

0.43(1) 

12(3) 

0.5(2) 

52 

-3.5 

3.76(4) 

0.17(1) 

19(4) 

0.7(2) 

53 

-5.0 

7.64(6) 

0.051(2) 

24(6) 

0.9(3) 

54 

-2.0 

0.91(4) 

0.49(1) 

9(4) 

0.3(2) 

55 

-3.5 

2.29(7) 

0.26(2) 

18(10) 

0.7(4) 


Acceptance v.s. <AH> 


6. CONCLUSIONS 
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Figure 1. Swap Acceptance Rates 


Ak = K 2 — Ki is increased. Figure shows the 
acceptance rate (A) as a function of (ATi). The 
dashed line is the graph of (^. It can be seen 
that the measurements and the predictions of the 
acceptance rate agree very well. 

We note that the simulations with parameters 
given by action matching technology (S'! — S'3) 
have a lower acceptance rate than the others. 
This is because the fluctuations in AH are larger 
using our pseudofermionic Hamiltonian, than in 
the action used to perform the matching. This 
issue is discussed more fully in Q. 

Column 5 in table § shows the integrated auto¬ 
correlation times of the plaquette for our simula¬ 
tions. Column 6 gives the corresponding ratios of 
PT to HMC autocorrelation times. These ratios 
are consistent with the predictions of the model 
described in section |[ 


We find that the acceptance rate of the PT al¬ 
gorithm drops very rapidly with Ak. This prob¬ 
lem is expected to get worse on larger lattices as 
AH, an extensive quantity, will have larger fluctu¬ 
ations. Over the range of available Ak the HMC 
autocorrelation times are equal within errors and 
the predictions of our autocorrelation model ap¬ 
ply. We estimate that connecting fast and slow 
decorrelating regions of parameter space would 
need a very large number of subensembles mak¬ 
ing PT impractical to use for lattice QCD with 
currently available computer technology. 
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